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Abstract. We revisit the method of cumulants for analysing dynamic light 
scattering data in particle sizing applications. Here the data, in the form of 
the time correlation function of scattered light, is written as a series involving 
the first few cumulants (or moments) of the distribution of particle diffusion 
constants. Frisken (2001 Applied Optics 40 4087, [T]) has pointed out that, 
despite greater computational complexity, a non-linear, iterative, analysis of the 
data has advantages over the linear least-squares analysis used originally. In 
order to explore further the potential and limitations of cumulant methods we 
analyse, by both linear and non-linear methods, computer-generated data with 
realistic ‘noise’, where the parameters of the distribution can be set explicitly. 
We find that, with modern computers, non-linear analysis is straightforward and 
robust. The mean and variance of the distribution of diffusion constants can be 
obtained quite accurately for distributions of width (standard deviation/mean) 
up to about 0.6, but there appears to be little prospect of obtaining meaningful 
higher moments. 
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1. Introduction 

Small particles of many kinds occnr natnrally in biological systems and the 
environment, and are nsed at some stage in many industrial processes |2]. 
Frequently the size of the particles is crucial to their function. Thus measuring 
particle size is important, and many techniques have been developed for this 
purpose |2]. One such technique, dynamic light scattering (DLS), can be applied to 
nanometric particles, size from a few nm to about 1 /im, which can be suspended in 
a liquid. DLS has the advantages of being quick and reproducible and of providing 
well-defined (though often limited) information about the particles. 

In DLS, coherent laser light scattered by a particle suspension forms a random 
diffraction pattern that fluctuates in time as the particles move in Brownian 
diffusion [SI a E]. For equal-sized particles, analysis of the time dependence 
of the scattered light yields the particles’ diffusion constant from which its 
size can be calculated. In the more common situation where the suspension 
is polydisperse — there is a distribution of particle sizes — DLS yields the 
Laplace transform of the distribution of diffusion coefficients. Thus, in principle, 
the latter, usually the quantity of interest, can be obtained by inverse Laplace 
transformation of the data. A variety of techniques have been used to perform 
this operation, including exponential sampling [6], regularization [7], maximum 
entropy [8], maximum likelihood [9] and non-negatively constrained least squares 
[To] . Because inverse Laplace transformation is particularly sensitive to inevitable 
experimental uncertainties in the data, these techniques are best suited to broad 
distributions (which may be multi-modal). Furthermore, the techniques can be 
quite complicated to operate and may require the input of prior information about 
the sample such as minimum and maximum particle size. A survey and critique 
of these methods, 20 years old but still valuable, was given by Finsy [TT] . 

A simpler approach to DLS data analysis, which was in fact the hrst one to 
go beyond fitting a single exponential, is the so-called method of cumulants m 
ESIE]. This provides a few lower moments, or cumulants, of the distribution of 
diffusion coefficients, and is the topic addressed in this paper. In the early days of 
DLS, when computing power was limited, cumulant analysis used a linear fitting 
method which did not require an iterative program. Later Frisken [1] pointed out 
that a much more versatile non-linear, iterative, fitting procedure is possible with 
modern computers. Frisken mm and others [121 demonstrated the value of this 
approach to analyze real experimental systems. Here we use realistic computer¬ 
generated data, where the parameters of the size distribution can be set explicitly. 
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to evaluate further the potential and limitations of non-linear cumulant analysis. 

In the next section we describe the background to the current situation. While 
much of this material has appeared before, we believe that it is helpful to provide 
a coherent account here. Section describes generation of the synthetic data and 
the analysis methods used. In section]^ we present the results which are discussed 
in section We hnd that, if the distribution of diffusion coefficients is not too 
broad, non-linear cumulant analysis offers a straightforward and robust method for 
determining its mean and variance. However, the prospects for obtaining higher 
moments are not promising. 


2. Background 


General references for this section are [31 m na El E]. DLS measures the 
normalised time correlation function g^'^^r) of the scattered light intensity J: 




(J)2 ^ 


( 1 ) 


where r is the correlation delay time. In many cases the intensity correlation 
function can be written in terms of the correlation function g'^^^r) of the scattered 
light held through the so-called Siegert relation 


= 1+P 


( 2 ) 


where (3 is the coherence factor, determined largely by the ratio of the detector 
area to the coherence area of the scattered light; /3 is usually regarded as an 
unknown parameter to be htted in the data analysis. In the simplest case of a 
dilute suspension of identical spherical particles in Brownian motion, g^^^r) is 
given by 

= exp(-rT) (3) 

where the decay rate V is 

r = Dq\ (4) 

D is the translational diffusion coefficient of the particles and q is the scattering 
vector (set by the scattering angle 9 and the wavelength A of the light in the sample 
through q = (dvr/A) sin(6'/2)). In turn, for spherical particles, D is given by the 
Stokes-Einstein relation 


QnrfR 


( 5 ) 
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where is the thermal energy, rj the viscosity of the solvent and R the particles’ 
radius. Equations ([T]-[^ form the basis of particle sizing by DLS. 

When a sample is polydisperse, containing particles of different sizes, each 
species gives rise to its own exponential decay in the held correlation function so 
that 

gW(r) = yG(r)exp(-rT)dr, (6) 

where G(r) is the normalized distribution of decay rates E (f G(r)dr = 1). Thus 
obtained from the measurement of through Eq. ([^, is the Laplace 

transform of G(r). In principle, therefore, the latter can be obtained by inverse 
Laplace transformation of the data. In practice, inverse Laplace transformation is 
an ill-conditioned problem in the sense that it converts small uncertainties in the 
data into large uncertainties in the recovered G(r) [18]. Put another way, unless 
there is a wide spread of particle size, the sum of exponentials implied by Eq. (§ 
looks not too different from a single, average, exponential. 

This limitation can be recognized and exploited by writing 


exr 


where F is the mean value of G(r), 
T = [ TG{T)dT, 


^(b(r) = exp(-rr) 


1 


where 


/i„ = / (r - frGiT)dT 


f)^], 

(7) 


(8) 

Eq. (g, 



(9) 


(10) 


are the central moments of the distribution of decay rates (the moments about 
the mean). Equation ([^ shows clearly how DLS data can be represented by an 
average exponential with correction terms that depend on G(r) (and hence on 
the particle size distribution). For a reasonably narrow distribution of decay rates 
and a reasonable range of scaled delay time f r, the higher-order terms in Eq. § 
become increasingly unimportant. Rewriting Eq. ([^ in terms of scaled time. 


(7(')(r) = exp(-rr) 


1 I ^ ^ h'S /p \3 I ^ 


( 11 ) 
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shows that /i2/r^, the normalized variance of G(r), is the simplest measnre of the 
departnre of g^^^r) from a single exponential. With Eq. ([^, Eq. (§ can also be 
written as 

In - 1 = ^ln/?-fr + - ..., (12) 

showing further how non-exponentiality appears as departure from linearity in a 
semi-logarithmic plot of the data (see Fig. |^. 

The original method of cumulants [I2l [131 [H] follows Eq. (12): the left-hand 
side, calculated from the data, is fitted to a polynomial of a few terms in delay time 
r, hence providing estimates of /3, T, etc. This method has the advantage that 
least-squares htting to a polynomial which is linear in the unknown coefficients is 
a soluble problem that does not require iteration in the computer program Hi. 
A disadvantage of the method is that, to keep the higher-order terms in Eq. (12) 
small, the data have to be truncated at around Tr = l(i.e. only data for Tr < 1 are 
kept) and it is not straightforward to determine the optimum truncation. There 
is, in fact, a trade-off between large random errors in the htted parameters if the 
data are truncated at too small a value of Tr and large systematic errors (but 
smaller random ones) if too much of the data are used m. 

Later, following rapid development of computer power, Frisken [T] pointed out 
that iterative, non-linear htting of 

2 


= 5 + ^ Jgxp(-rr) 


1 + 


(13) 


obtained from Eqs. (|^ and 0 . is a more robust procedure. This approach has 
several advantages over the linear method. First, it is not necessary to truncate 
the data since the divergence of the higher-order terms in Eq. (13) is suppressed 
by the decaying exponential pre-factor. Second, the method allows the “baseline” 
B to be regarded as a parameter to be htted. In an ideal experiment, B should 
be 1 (Eq. (|^). In practice B can diher slightly from 1. For example, slow drift 
of the laser intensity or of the gain of the detector leads to a spurious correlation, 
B > 1, which is almost independent of time over the span of the data. 

Frisken [T] used Eqs. (12) and (13) to analyze experimental data and 
clearly demonstrated the advantages of the non-linear method outlined above. 
Subsequently Hassan and Kulshreshtha (TB] performed a similar analysis of 
experimental data and also considered simulated data for known distributions of 
decay rate G(r). However they only included terms up to second order in time 
and their simulated data did not take account of the uncertainty (noise) that is 
inevitable in an experiment. In this paper we compare the two methods of analysis 
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using simulated data with realistic noise over a range of polydispersities, with the 
aim of determining more precisely the potential and limitations of the non-linear 
method. 

A note on terminology: Koppel na pointed out that, formally, the logarithm 
of Eq. @, is the cumulant generating function [20] for the distribution 


G(r). Thus the expansion of this quantity, Eq. (12), is a power series in which 


the coefficients are the cumulants of G(r); it is from this observation that the 


linear analysis based on Eq. (12) gets its commonly-used name, the “method of 
cumulants”. However, in the non-linear analysis using Eq. (13) emphasized here. 


it is the central moments rather than the cumulants that are relevant; thus what, 
to follow the custom, we have here called non-linear cumulant analysis might more 
logically be called the “method of moments”. (In fact, comparing Eqs. (12) and 


(13), we see that the cumulants only differ from the central moments at order 4 


and higher rao 


3. Methods 

Mainly for mathematical convenience, we assume a Schulz distribution of decay 
rates: 




(14) 


this is a two-parameter distribution defined by mean decay rate T and width a 
(normalized standard deviation), given by 
p_f2 1 




(15) 


n z +1 

Sample plots of the Schulz distribution are shown in Fig. ([^. 

Substitution of Eq. ( |I4) ) into Eq. (|^ gives 

= (l + n^fr) ; ( 16 ) 

(series expansion verifies that Eq. ( p)6| reduces to a single exponential, Eq. ([^, 
as the width a of the distribution tends to zero). The moments about the origin 
of the Schulz distribution are 


T" = y T^G{TdT) = r”(l + (n - l)a2)(l + (n - 2)a^) ... (1 + cr^), (17) 

giving, for the central moments, Eq. (0. 

^■2 2 r, 4 j ^-4 Q /' 4 I o 6\ 

pj = CT , P = 2ct and pj = 3 (ct + 2ct ) ; 


( 18 ) 
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Figure 1. The Schulz distribution, Eq. (14), for indicated values of standard 


deviation a. 


the fourth cumulant of the distribution is 


^■4 „ 

f4 




(19) 


Synthetic “data” g^'^^r) were constructed for a range of distribution widths 


0 < cr < 1 from 


9 


= 


B + /3[s™(t)]' 


( 20 ) 


with g^^^r) given by Eq. (16); here B is the baseline, cf. Eq. (13). To mimic 
modern photon correlators, 150 data points were logarithmically spaced in scaled 
delay time Er in the range 10“^ < Er < 10^. To mimic experimental noise, 
each value of was multiplied by a random number drawn from a Gaussian 

distribution of mean 1 and standard deviation s. For most of the analysis we took 
s = 10“^, corresponding to an uncertainty of one part in a thousand on each data 
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point. This is the typical magnitude of counting errors in a DLS experiment. We 
also looked briefly at noisier data, up to s = 10“^. For each value of a, 20 data 
sets with different random noise were analyzed, allowing calculation of the means 
and standard deviations of the htted parameters. 

Two analyses of the data were performed. In the standard cumulant analysis, 
we assumed B to take its ideal value of 1 in Eq. (20). Then In \/— 1, 
calculated from Eq. ([^, was htted by linear least squares to Eq. (12) [19]. Fits 
to polynomials in r of order one, two, three and four were performed, providing 
estimates of an increasing number of the moments /r„. For simplicity, both F and 
(3 were assumed to be 1 when generating the data, but were taken as parameters to 
be htted in the analysis. In this linear cumulant analysis, the data were truncated 
when — 1 dropped to 10% of its initial value. 

In the second, non-linear, analysis, the simulated data for were htted to 

Eq. (13) using a variable metric method As with the standard analysis, four 
orders of ht were performed. The data were not truncated and the background B 
was regarded as an additional boating parameter. 


4. Results 

Figure a) shows an example of data htted successfully by the fourth-order non¬ 
linear procedure. Input values were a = 0.4 and B = 1.01, and all the htted 
parameters are within the expected uncertainty of the input values. The dashed 
line shows a single exponential with decay rate F = 1. The hgure illustrates 
how even a signihcant spread of particle size (see Fig. [^, cr = 0.4, leads to a 
correlation function that does not diher much from a single exponential. In the 
semi-logarithmic representation of Fig. [^b), the diherence, at larger delay times, 
is more apparent. 

4.1. Mean decay rate 

Figure 1^ shows the deviation of the htted mean decay rate F from its input value, 
1, as a function of polydispersity a for both linear and non-linear hrst- to fourth- 
order hts. A ht of order 1 is the equivalent of force-htting the data to a single 
exponential. It is clear that, as soon as polydispersity becomes signihcant, the 
hrst-order hts seriously underestimate F. However, adding just one parameter, 
/X 2 , in the second-order hts immediately allows reliable estimates of f up to 
polydispersities of 0.4 to 0.5; fourth-order hts estimate F reliably over almost 
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Figure 2. Example of a fit to simulated data, (a, top) Crosses: simulated 
data, Eqs. (20) and (16), with B = 1.01,/3 = l,f = l,cr = 0.4 and noise 
s = 10“^. Solid line: fourth-order non-linear fit of the data (Eq. 
Residuals, difference between data and fit, are indicated. Dashed line; single¬ 
exponential decay with E = 1. (b, bottom) Semi-logarithmic representation of 
the same data. 
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Figure 3. Deviation of the fitted mean decay rate F from its input value, 1, 
as a function of polydispersity a for both linear and non-linear first- to fourth- 
order fits. The error bar indicates the standard deviation of the fitted f for the 
fourth-order non-linear fit across the 20 generated data sets at cr = 0.4. See text 
for discussion. 


the whole range of polydispersity considered. The error bar in Fig. [^indicates the 
standard deviation of the htted F for the fourth-order non-linear £t across the 20 
generated data sets at a = 0.4, and shows that F can be obtained with an accuracy 
of better than 0.5% for moderately polydisperse samples. 

In the estimation of F for moderately polydisperse systems, there is little to 
choose between non-linear and (truncated) linear hts. Intriguingly, at third-order, 
the non-linear £t actually does worse than the linear one. We note that, for a 
symmetrical distribution of decay rates, the odd-order central moments are zero. 
Thus, in general, /i 4 r^ can be larger than in Eq. (13) even at Fr < 1, 


and 


there is no justihcation for a £t that includes pa but not / 14 . 
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Figure 4. Deviation of the fitted second central moment y ,2 from its input 
value (T^ (Eq. (18) with F = 1) for linear and non-linear fits up to fourth order. 
Similar to Fig. the error bar indicates the standard deviation of the fitted ^2 
for the fourth-order non-linear fit across the 20 generated data sets at cr = 0.4. 


Under the same conditions as for Fig. i Fig. a shows the deviation of the 
htted second central moment /i 2 from its input value (Eq. (18) with F = 1). 
For this parameter, all three non-linear hts do much better than the linear ones, 
giving accurate estimates of /i 2 for polydispersities up to 0.3. The fourth-order 
non-linear £t gives a good estimate of /t 2 up to about a = 0.6. At a = 0.4, where 
/i 2 = = 0.16, uncertainty in the determination of fi 2 is about ±0.02, a hnding 

consistent with experience in experiments [H] . These values translate into a quite 
accurate determination of the width, cr(£tted) = 0.40 ± 0.025 
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Figure 5. Deviation of the fitted third central moment from its input value 
2 ( 7 ^ (Eq. ra with f = 1). 


Figure shows the deviation of the third moment from its input value 2(T^ 
(Eq. (18) with F = 1). We have mentioned above that there is no justihcation 
for a third-order £t. At hrst sight the fourth-order £t appears to do well up to 
a = 0.4. However we note that, for the fourth-order fit at a = 0.4, the uncertainty, 
±0.045, in the recovered pa is barely smaller than the value, 0.0512, of ps itself. 
Thus, at least for this distribution and noise level, DLS can do little more than 
hint at the sign of ps. 


4-4- Fourth moment 

We hnd that, up to about a = 0.4, the uncertainty in the htted value of /i 4 is 
comparable to its actual value. On the other hand, above a = 0.4, as /i 4 becomes 
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larger, its fitted value is consistently lower than its actual value. Thus, as with /is, 
we obtain little useful information. 


4-5. Baseline 

Provided that the data extend to long enough times, as in Fig. we find that the 
baseline B is recovered accurately by the non-linear hts. Furthermore, treating 
the baseline as a parameter to be htted introduces almost no further uncertainty 
in the other fitted parameters (compared to the case where the baseline is fixed at 
its input value). 


4-6. Initial guesses 


The non-linear htting procedure, Eq. (13), appears to be very stable with respect 
to the choice of initial values for the htted parameters. The simplest approach 
is to obtain initial values of fl and F from a hrst-order linear £t of a short-time 
portion of the data, to assume that the initial value of the baseline i? is 1, and to 
take initial values for the higher moments, /i 2 ,/^ 3 , etc., to be zero. 


4-1 ■ Larger noise 

With noisier data, s = 10“^ (see section 3), fourth-order hts were not useful, giving 
large uncertainties in both F (several percent) and fi 2 (typically ±0.20). However 
second-order non-linear hts gave F to within about 2% with an uncertainty in /i 2 
of about ±0.05, meaning that polydispersities a (= y^) greater than 0.20 to 0.25 
could still be detected. 


5. Discussion 


Several conclusions may be drawn from these results. First, we hnd that non¬ 
linear hts are more accurate and more straightforward to perform than linear hts 
(where a somewhat arbitrary truncation of the data is necessary). This is seen 
clearly in Fig. which shows that non-linear hts return a much smaller systematic 
error in the second moment /i 2 than linear hts. Second, we have pointed out that 
there is no justihcation for performing third-order hts, either linear or non-linear: 
the r^-term in Eqs. (12) and (13) can, in general, be comparable to or larger 
than the r^-term even at small Fr. Third, performing a fourth-order, rather 
than a second-order, ht reduces the systematic error in the htted /i 2 (Fig. |^. 
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However, even for polydispersities as large as a = 0,6, fourth-order fits do not 
provide reliable estimates of the higher moments, and /i 4 . This hnding agrees 
with experience: we are not aware of any experiment where /is or /i 4 have been 
convincingly measured. 

It might at hrst appear surprising that including unknown parameters /is and 
/i 4 in the non-linear £t should improve the htting of /i 2 even though the determined 
values of /is and /i 4 are themselves unreliable. Comparisons of the surface |I9] 
in response to combined changes of /i 2 with f and /i 2 with /is reveal that the 
additional £t parameters signihcantly reduce the correlations between /i 2 and T 
i.e. the minimum area in the surface becomes perpendicular to either the /i 2 
or the r axes when /is and /i 4 are included in the £t. Without these parameters 
the minimum is typically aligned parallel to /i 2 = T making it more difficult to 
determine the optimal £t values. 

Despite the somewhat negative conclusions of the first paragraph of this 
section, we are able to suggest a robust and relatively straightforward procedure 
to obtain the hrst two moments of the distribution of diffusion constants, T and 

p2'- 


111 


IV 


Set up a photon correlator with channels logarithmically-spaced in time so 
that the measured intensity correlation function extends well into the 

long-time baseline B] run the experiment for long enough that the typical 
noise on a data point is about 1 part in 10^. 

Obtain initial estimates of the coherence factor f3 and mean decay rate T from 
a short-time linear £t of the data, and take other initial estimates as baseline 
5 = 1, higher moments P 3 , Pi = 0. 

Perform a fourth-order non-linear £t of g^‘^\T) to Eq. (13). If the htted 
second moment /i 2 /f^ is less than about 0.40 (corresponding to a ~ 0.6 in 
Fig. 4), then the htted f and /i 2 /f^ should be good estimates of the true 
values, accurate to around 0.5% and ±0.02 respectively. 

If, in (hi), the htted second moment /i 2 /f^ is less than ~ 0.1 (polydispersity 
a < 0.3), then more precise estimates of the values of P and /i 2 /r^ should be 


obtainable from a second-order non-linear ht of g^'^^r) to Eq. (13) (see Figs. 
3 and 4). 

Here we have only considered a Schulz distribution of decay rates G(r) (Eq. 
(14)). However it is straightforward to show that the moments of a distribution 
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can be written [ 221123 ] 


pn 


1 + 


n(n 


1 ) 


-a 


+ Oia^) 


( 21 ) 


For, fairly narrow, fairly symmetrical distributions it is sufficient to keep only the 
first two terms so that the higher moments can, to a good approximation, be 
written purely in terms of the standard deviation a. Thus we may expect that the 
general discussion above will apply to any distribution which is narrow enough. 
For the Schulz distribution, we have shown here that “narrow enough” means a 
less than about 0.6; but, for more skewed or flatter (larger kurtosis) distributions 
than the Schulz, this upper limit might need to be reduced somewhat. Repeating 
the programme of this paper for other distributions such as the lognormal or a two- 
component mixture would provide more quantitative conclusions, but the general 
picture is unlikely to change. 

The method of cumulants described here is relatively simple to implement and 
does not require the input of any prior information about the sample. It should 
therefore be the first approach of any experimenter presented with an unknown 
sample. When this hrst measurement returns a large enough second moment, it is 
probably worth trying some of the more complex analysis methods listed in section 
1. Experience with these methods HD shows that they have difficulty resolving 
bimodal distributions when the ratio of the two sizes is less than about 3. This 
corresponds to a standard deviation a ~ 0.5 and a second moment ~ 0.25. 

In general, though, these methods work best for signihcantly broader distributions 
than this. 

In this paper we have limited consideration to analysis of a single DLS 
measurement. For completeness, we mention that, for particles large enough that 
there is significant angular variation in the intensity of light that they scatter 
(radius greater than about 50 nm), performing a combined analysis of data 
taken at different scattering angles can frequently provide much more detailed 
information. For example, by measuring the angular dependence of the apparent 
average diffusion constant of particles which show a minimum in their angular 
intensity profile, it is possible to measure very small polydispersities [23] • R 
has recently been demonstrated that a Bayesian analysis of the full correlation 
functions measured at several scattering angles can be very powerful in resolving 
multi-modal distributions [25] . 

We note that obtaining information about the distribution G(F) of decay 
rates (or diffusion constants) is usually not the ultimate goal of an analysis of DLS 
data; rather, one is interested in the distribution of particle sizes. Because the 
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contribution of each particle species to G(r) is weighted by the intensity of light 
scattered by that species and because the particle radius is, through Eqs. Q and 
([^, inversely proportional to the decay rate F, obtaining the size distribution from 
G(r) is not straightforward. Nevertheless, for many systems, such as solid spheres 
|1] , spherical shell-like particles |T5] and random-coil polymers |ll] , the measured 
moments F and /i 2 /r^ can be related to moments of the particle size distribution. 
For example, for homogeneous spheres much smaller than the wavelength of light, 
it can be shown 01121] that the effective particle radius obtained by substituting 
the measured f in Eqs. (HI and ([^ is 

ksTg^' 


R. 


eff 




= 


and that 


fi2_ R^R^ ^ 

F2 (7^)2 


( 22 ) 


(23) 


where the R"^ are moments of the distribution of particle radii. For narrow size 
distributions, where an approximation like that of Eq. (21) can be applied, these 
results reduce to the simpler (and useful) expressions = i? (1 -|- 5cr|.) and 
/i 2 /F^ = where ctr is the standard deviation of the radius. 

Finally we comment that the conclusions reached here — useful determination 
of F and /i 2 , but little prospect of obtaining higher moments — are not too 
different from those reached by Koppel [T2| when introducing the linear cumulant 
method more than 40 years ago. However we have shown that non-linear htting 
gives a simpler and more robust procedure which avoids the necessity of a rather 
arbitrary truncation of the data; we have estimated the upper limit of reliability 
of cumulants methods at about /i 2 /F^ = 0.4, polydispersity a ~ 0.6; and we 
have suggested that if an initial analysis by the cumulants method yields a second 
moment more than about /i 2 /F^ = 0.25, it is probably worth trying to obtain more 
information from one of the more complex analysis methods, based on inverse 
Laplace transformation, mentioned above. 
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